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ABSTRACT 

There is an observational correlation between astrophysical shocks and non-thermal 
particle distributions extending to high energies. As a first step toward investigating the 
possible feedback of these particles on the shock at the microscopic level, we perform 
particle-in-cell (PIC) simulations of a simplified environment consisting of uniform, in- 
terpenetrating plasmas, both with and without an additional population of cosmic rays. 
We vary the relative density of the counterstreaming plasmas, the strength of a homo- 
geneous parallel magnetic field, and the energy density in cosmic rays. We compare the 
early development of the unstable spectrum for selected configurations without cosmic 
rays to the growth rates predicted from linear theory, for assurance that the system is 
well represented by the PIC technique. Within the parameter space explored, we do 
not detect an unambiguous signature of any cosmic-ray-induced effects on the micro- 
scopic instabilities that govern the formation of a shock. We demonstrate that an overly 
coarse distribution of energetic particles can artificially alter the statistical noise that 
produces the perturbative seeds of instabilities, and that such effects can be mitigated 
by increasing the density of computational particles. 

Subject headings: cosmic rays, instabilities, plasmas, shock waves 
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Introduction 



Shocks form in a wide variety of astrophysical environments, from planetary bow shocks in 
the heliosphere to colliding clusters of galaxies. The presence of nonthermal partic le populations i n 
some of these environments is inferred from radio, X-ray, and 7-ray observations ([Reynolds! 120081 ) . 
and a number of theoretical mechanisms link stro ng shocks (such as those found in young supernova 
remnants) to particle acceleration processes fe.g.. lDrurvill98.4 IWebb et al.lll983l ) . An understanding 
of the nonlinear coupling between the shock, the energetic particles, and the spectrum of any excited 
waves is essential to the proper interpretation of the observational signatures of shocked systems 
as well as to a correct understanding of the origin of cosmic rays. 

In most cases the density of the shocked medium is sufficiently low for collisions between par- 
ticles to be infrequent; such collisionless shocks are mediated through collective electromagnetic 
effects. The thickness of the shock transition region and the nature of the microphysics occurring 
therein may play a key role in the injection efficiency of t hermal part icl es into the accelera tion 
processes thought to accelerate cosmic rays to high energy iiellll978J lbl: IScholer et~aH hood ) . In 
addition, efficient particle acceleration associated with some shock environments raises the question 
of the extent to which a significant cosmic-ray contribution to the local energy density can influence 
the structure of the shock front itself. This may be of particular importance in starburst galaxies, 
where massive stars and frequent supernovae increase the availability of sites for particle acceler- 
ation. The emission of TeV 7-rays from some nearby starburst g alaxies suggests an abundance of 
cosmic rays that significantly exceeds the values observed locally (jAcciari et al.l l2009). 



A number of mechanisms exist whereby cosmic rays may play a role in shaping the shock and 
its environment. A significant cosmic-ray contribution to the pressure upstream of a shock is un- 
derstood to result in substantial modification to the shock environment (|Eichlerlll984l ); this effect is 
expected to produce an effective shock compression ratio that is perceived as bei ng larger by cosmic 
rays of higher energy, possibly hardening the local cosmic-ray source spectrum (jMalkov and Drurv 
200ll ) . Under some conditions, the current carried by cosmic-ray ions diffusing ahead of the shock 
may excite instabilities in the upstream interstellar medium that lead to large-amplitude magnetic 



turbulence (Bell 2004: iLemoine and Pelletier 



2011 



Rabinak et al 



2011 



Nakar et al.ll201ll ) or other 



substantial alteration s to the upstream environment, such as the heating of thermal electrons (e.g., 
Rakowski et al.ll2008l ). However, these effects occur on spatial scales much larger than the thickness 
of the subshock and may operate quite independently of the processes occurring therein. 

In this paper, we turn our attention to those processes within the subshock. In order to 
quantify or constrain the effect of a "spectator" population of cosmic rays (that is, any cosmic rays 
present that have already been accelerated, locally or elsewhere) on the formation and evolution 
of nonrelativistic collisionless subshocks, we design simulations to explore the initial linear and 
subsequent nonlinear growth of instabilities that contribute to shaping the shock transition region. 
We restrict our focus to the interpenetration layer, in which two counterstreaming plasma shells 
overlap in space. This is fertile ground for the growth of well-known and well-studied instabilities 
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(jSilva et al.l 120031 ; iFrederiksen et al.l 12004 ; iBretl 120091 ). and the two-plasma system will respond 



accordingly. We then repeat our simulations with an added cosmic-ray component, a plasma 
consisting of highly relativistic particles, so we may compare the systems at both the early and late 
stages of their evolution and identify whether the presence of the cosmic rays has any effect. 

The scale of the instabilities dictates the computational approach to modeling them. Hydro- 
dynamical models of shocks are incapable of accurately resolving length scales much smaller than 
a particle mean free path and must therefore approximate subshocks as discontinuities, which is 
inadequate for our purposes. A self-consistent kinetic approach is necessary for a more complete 
exploration of the relevant small-scale physical processes, but simulating the entire subshock thick- 
ness is prohibitively costly. We elect a simulation that represents only a homogeneous portion 
of the subshock interior, and therefore we cannot account for a number of instabilities that arise 
from spatial gradients or particle reflection. Here we use particle-in-cell (PIC) simulations in two 
spatial dimensions to model the interaction of counterstreaming plasma flows in the presence of a 
third, hot plasma of energetic particles. The PIC technique is a (particle-based) kinetic approach 
to modeling the self-consistent evolution of an arbitrary distribution of charged particles and elec- 
tromagnetic fields and thus is well suited to the nonlinear development of unstable plasma systems. 
In order to isolate the subshock instabilities from larger-scale spatial effects such as those arising 
from systematic charge separation, all distribution functions are initially homogeneous with respect 
to position. We consider the simplified case in which two electron-ion plasmas, not necessarily of 
equal density, move nonrelativistically through each other. This movement may also be parallel to 
a uniform magnetic field. We first observe the evolution of the drift velocities, field amplitudes, 
particle distributions, and wave spectra in the case when cosmic rays are negligible. We then in- 
clude cosmic rays whose energy density is an order of magnitude below the kinetic energy in the 
bulk flow of the plasmas. We also explore the case of an artificially high energy density in cosmic 
rays, exceeding the value of equipartition with the bulk flows, to gain insight into effects that may 
be too small to arise in the case intended to represent a more realistic environment but that may 
nevertheless appear in regions where energetic particles are unusually abundant. However, our 
results suggest that even a high abundance of cosmic rays is not sufficient on its own to produce a 
significant deviation from the usual evolution of the instabilities shaping a shock. 



2. Objectives and approach 

To constrain the extent to which cosmic rays might influence the physics of the shock transition 
layer, we perform a series of simulations. Using the benchmarks of system behavior outlined in 
Section 12.31 we characterize the response of the counterstreaming-plasma system to the cascade 
of instabilities arising from interactions among the various plasma components in the absence of 
cosmic rays. We then repeat the simulations with cosmic rays present, so as to facilitate the 
side-by-side comparison of the various benchmarks. The parameter choices for the systems that 
we explore in this manner are described in detail in Section 12.11 and our simulation model and 
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its implementation are described in Section 12,21 Then, a test of selected parameter configurations 
against the predictions from analytical studies of related systems is described in Section Following 
a brief discussion, we conclude in Section [5j 

2.1. Parameter-space configurations 

The physical system under consideration is modeled as one homogeneous electron-ion plasma 
moving relative to another. These plasmas may in general have different densities. In addition, 
a uniform magnetic field may be aligned with the flow direction. Finally, a population of cosmic 
rays may be present, at rest in bulk in the center-of-momentum frame of the two plasmas. This is 
the reference frame chosen for the simulation, so the two plasmas flow in opposite directions; we 
adopt the nomenclature of "stream" and "counterstream" to distinguish between them. By our 
convention, the stream's velocity is in the — x-direction, antiparallel to the guiding magnetic field 
when one is present. 

Our primary simulations are sensitive to variations in three parameters: the density ratio 
between the stream and the counterstream, the strength of the guiding magnetic field, and the 
energy density in the cosmic rays. 

The inter-plasma density ratio w = n s /n cs takes three values and their respective designations: 
the "symmetric" case w = 1, the "intermediate" case w = 0.3, and the "dilute" case w = 0.1. The 
velocity of the stream is fixed at v s = — 0.2cx, while the velocity of the counterstream obeys the 
relation n cs v cs + n s v s = 0; in the dilute-stream case, the relative flow speed is therefore 0.22c. 
Likewise, the density of the counterstream is fixed, so the stream density alone varies; the total 
electron density n e and thus the electron plasma frequency w pc = y e 2 n e /eom e (where eo is the 
vacuum permittivity, and n e = n s +n cs is the cumulative electron density) is largest in the symmetric 
case and reduced by a factor yTT/2 in the dilute case. 

The magnetic field Bq )X may be absent, or present at either of two amplitudes given by the 
ratio of the electron cyclotron frequency f2 e = eBo tX /m e to the electron plasma frequency, b = 
Q e /ujp e . The "absent" magnetic field refers to b = Bq <x = 0. The values designated "weak" and 
"strong" correspond to b = 0.01 and b = 0.1, respectively; in a plasma of electron density of, 
e.g., n e ~ 1 cm~ 3 , a magnetic field of ~ 300 fiG is necessary for b = 0.1. Since the speed of 
Alfven hydromagnetic waves va = be/ Wl + mi/m e is at most 0.014c for our choice of m,i/m e = 
50, all of the plasma collisions we consider have Alfvenic Mach numbers significantly larger than 
unity (up to oo in the case when 6 = 0). Note that because the electron plasma frequency is 
not independent of the density ratio w in our simulations, neither is the absolute magnetic-field 
amplitude corresponding to a particular value of b. 

All simulations include cosmic-ray particles consisting of electrons and ions, initialized ac- 
cording to an isotropic distribution function and a single speed (Lorentz factor 7cr = 50) whose 
statistical weight wcr = ncR/n s is adjusted to three levels: "negligible" when w;cr7cr = 10 -8 , 
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"present" when u>cr7cr = 10 -3 , and "abundant" when u>cr7cr = 10. Since wqr is defined 
in terms of the stream density n s , the absolute energy density in cosmic rays also varies with 
the density ratio w, being 10 times larger in the symmetric case than the dilute case for each 
value of u>cr- Neglecting the contribution of electrons, the bulk kinetic energy in the stream 
plus counterstream is of order n cs miVgW(l + w)/2, while the cosmic-ray energy density is of order 
nCRTCR^jC 2 = 25wcR"fcRwn cs 'miVg , where we have used the relation v s /c = 0.2. Thus, the ratio of 
cosmic-ray energy density to bulk kinetic energy density is 50u?cr7cr/(1 + w): considerably larger 
than unity for the "abundant" case and a few percent in the "present" case. 



2.2. Simulation setup 



O ur simulations e mploy a modified version of the relativistic electromagnetic PIC code TRIS- 
TAN (jBunemanl Il993l ). updated for parallel use with MPI, operating in two spatial dimensions, 
with three-component velocity and field vector s (2D3V), with peri odic boundary conditions. The 
charge-conserving current deposit ion routine of lUmeda et al.l (|2003l ) and the field update algorithm 
with fourth-order accuracy from iGreenwood et al.1 (120041 ) are the most prominent additions, as 
well as digital filtering of electric currents to suppress small-scale noise via an iterative smoothing 
algorithm. 

The primary set of simulations, 27 in total, was conducted on a spatial grid in the x—y plane 
of size 280A se x 180A se (periodic boundary conditions in x and y, with elongation in the flow 
direction x), where A sc = c/uj pc = 10 A is the electron skin depth, set to 10 grid cells of length A. 
The electron plasma frequency w pe is determined from the sum of the stream and counterstream 
electron densities only and thus depends on w but not on wqr. Supplementary high-resolution 
simulations in which A se = 30A were performed on a grid of more cells but representing a smaller 

- 1 » 225t for the A se = 10 A 
665t for the high-resolution A se = 30A simulations. 



physical region, 128A se x 96A se . The time step St was chosen such that w pe 



simulations, or 



Six separate particle populations from three plasmas are modeled: the "stream" moving in 
the — x-direction, the "counterstream" moving in the +x-direction, and the energetic particles 
representing cosmic rays. Within each plasma, ions and electrons of charge ±e and mass ratio 
irii/m e = 50 have equal charge density and a common drift velocity so that the entire setup has 
no net current and no charge imbalance. The artifi cially low mass ratio expedites the simulations, 
but may change some modes (IHellinger et al.1 120071 ). Quite a few of these possibly modified modes 
are not included here, partly because we do not consider the spatial structure of subshocks, partly 
because they do not fit onto the computational grid, and partly because we do not simulate situa- 
tions involving an oblique or a perpendicular large-scale magnetic field. In any case, the analytical 
treatment of instabilities presented in Section [3] accounts for the small mass ratio, and in that sense 
our approach is consistent, at least for the linear phase. 



Each cell in a primary (high-resolution supplementary) simulation is initialized with a total 
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of 90 (120) computational particles: 20 stream ions, 20 counterstream ions, and 5 (20) cosmic-ray 
ions; and an electron for each ion. The physical density of each plasma is manipulated through the 
assignment of the appropriate statistical weights, w and u>cr, to the various particle species. 

The stream and counterstream, viewed from their respective rest frames at v s = —0.2cx and 
v cs = w x 0.2cx, are described by a Maxwell-Boltzmann distribution in which the electrons' most 
probable speed is given by v^ e = 0.01c and the ions are in equilibrium with the electrons. The 
cosmic rays, whose rest frame is the simulation frame vqr = 0, are isotropic and each is initialized 
with Lorentz factor 7cr = 50, regardless of whether it is an ion or an electron. Placing the 
cosmic-ray population at rest in the center-of-momentum frame mini mizes streaming in t he collision 



zone, and therefo re reduces known cosmic-ray-driven insta bilities ( Achterberj Il983 ; Bell 



Revil 



2009 



e et al 



Ohira et al 



2006) , wh i ch can be independently studied (e.g. 



2009; 



Stroman et al 



2009 



Niemiec et al 



Riquelme and Spitkovsky 



2008 



200 



2004 



Gargate et al 



Luo and Melrose 



We are interested in whether or not cosmic rays can modify instabilities operating at subshocks, 
and therefore suppress cosmic-ray streaming instabilities . 



2.3. Behavioral benchmarks 

To provide a basis for comparison among different cosmic-ray densities wqr for each combi- 
nation of plasma density ratio w and magnetic-field amplitude b, we select the following attributes 
of the system for study: the drift velocity of each particle population, the instantaneous root- 
mean-square amplitude of the parallel and perpendicular components of the electric and magnetic 
field over the entire simulation domain, the effective temperature of each stream or counterstream 
particle species, and the spectrum of excited wave modes in the magnetic field. 

As there is no initial bulk motion in the perpendicular directions, only the parallel component 
V x of drift velocity is considered. The electric field amplitudes are presented in units of the scaling 
electric field E u = u) pe cm e /e (equivalent to the field at which the electric energy density is half 
the electrons' rest-mass energy density, eoE 2 /2 = N e m e c 2 /2); the magnetic field multiplied by c is 
expressible in the same units. 

The particle distributions may not remain strictly Maxwellian throughout the duration of the 
simulation. As a surrogate for temperature, therefore, the mean random kinetic energy of the 
electrons and ions of the stream and counterstream is calculated by determining the systematic 
velocity component within a 10A x 10A region and eliminating this local bulk motion via an 
appropriate Lorentz transformation; the mean post-transformation Lorentz factor 7' corresponds 
only to the random motion. 

Finally, we will explore the effect of cosmic rays on the time evolution of the spectral decom- 
position of the perpendicular (out-of-plane) magnetic field B z into its spatial Fourier components, 
both parallel (wave number k x = kn) and perpendicular to the drift (wave number k y = k±). 
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3. Comparison with analytical beam plasma predictions 

As a test that our simulation results were consistent with theory, we applied the methods 



of iBretl (|2009l ) to selected stream-counterstream configurations without cosmic rays. Whether 
magnetized or not, beam-plasma systems (in which a fast, dilute "beam" plays a role comparable 
to that of our stream, with the dense "plasma" representing our counterstream) are susceptible 
to a host of both electrostatic and electromagnetic instabilities. For flow-aligned wave vectors, 
electrostatic modes such as two-stream or Buneman are likely to grow. In the direction normal 
to the flow, the filamentation instability (sometimes referred to as "Weibel") is usually excited as 
well. Finally, modes with wave vectors oriented obliquely are likewise unstable, so that the unstable 
spectrum is eventually at least two dimensional. 

The full spectrum has been first evaluated solving the exact dispersion equation in the cold 
approximation, accounting thus for a guiding magnetic field as well as finite-mass ions. It turns 
out that for the present configuration, ions play a very limited role (in the linear phase) and are 
not responsible for any unstable modes which would not be excited if they were infinitely m assive. 



Unlike settings exhibiting nonresonant modes, for example ([Bell 2004: IStroman et al.ll2009l ). where 
a single proton beam is considered without electrons moving at the same speed, we are here dealing 
with a plasma-shell collision, where protons and electrons are comoving. As a result, the effect 
of finite-mass ions is simply a first-order correction to the electronic spectrum. The dispersion 
equation displays the very same branches, and the growth rate is altered by a quantity proportional 
to 0(m e /rrii). 

The dispersi on equation for arbitrarily oriented modes with the flow along the x-axis reads 
(|Bret et alJbnid ) 

{u?e xx - kyC 2 )(uj 2 e y y - k 2 x c 2 ) - {u?e xy - k y k x c 2 ) 2 = 0, (1) 



where the tensor elements e a p are given by 



^ uj 2 J 7(P)%i j w 2 J 7(PJ mju - k • p/7(p) 

and the sum runs over the species involved in the system. For each species, the distribution function 
fj includes both the stream and the counterstream, and the corresponding plasma frequency co p j 
is calculated using the sum of their densities. The results in the cold-plasma limit are evaluated 
considering Dirac's delta distribution functions. Though lengthy, calculations are straightforward. 
Setting k y = kj_ = in the equations above allows to derive the dispersion equation for flow- 
aligned modes such as the two-stream ones. Setting k x = kn = gives the dispersion equation for 
the filamentation modes. We present analytical results for the symmetric and the diluted beam 
cases. While such expressions can be derived considering either w = 1 or w -C 1, expressions valid 
for any density ratio w are much more involved (when they exist). This explains why the results 
below for w = 1 cannot be derived from those with w <C 1. 



2 ,. FtfO , ,2 



k-(M 
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For flow-aligned wave vectors, the most unstable wave vector k\\ m and maximum growth rate 
7 m read (with evaluation corresponding to the configuration displayed in Figure [I]) 



h x Av + m e /m, _ 

Vil s 



C 



7m y/l + me/mi 

— — (J. 68, symmetric case, (oJ 



a;- 



pc 



V2T 



3/2 



and 



Aw 

h m X se — = vl + w = 1.04, it; < 1 



c 



7m V3(l + w) ^(l + me/m,) 1 / 3 

= o4/3 p = °- 31 > dllute case > ( 4 ) 

where r s = ( 1 — t> s 2 /c 2 ) is the bulk Lorentz factor of the stream, which moves in the simulation 
frame with speed v s = 0.417c/(l + w). As seen in Figure [TJ panels (a) and (f), oblique modes 
dominate for the mildly relativistic conditions of the simulation. The full-spectrum maximum 
growth rate is thus slightly larger than the numerical values calculated above for modes propagating 
along the flow. 

For wave vectors normal to the flow, the growth rate reaches its maximum for k± = oo, with 




7m v s 1 + m e /mi 

-> — — = 0.41, symmetric case, (5) 



and 



7m v s w(l + w)(l + m e /mi) 

= — \ = 0.12, dilute case. (6) 

w pe c y T s 

These filamentation data have been calculated neglecting the magnetic field, which is small (f2 e = 
O.Olcjpc) in the present setup. Electrostatic unstable modes propagating along the flow are rigor- 
ously insensitive to the flow-aligned magnetic field. As long as m e /mj <C 1, the two-dimensional 
linear spectra computed with or without finite-mass protons are indistinguishable. The hot spectra, 
accounting for the Vth,e = 0.01c thermal spread for the electrons, have then be en calculated con 



sidering infinitely heavy protons and using the fluid approximation described in iBret and Deutsch 



(2006). 



The predicted growth rates for the two-dimensional k^,k± wave- vector space are plotted in 
Figure [lj panels (b) and (g) for the dilute and symmetric streams, respectively, when the stream 
velocity is 0.417c (corresponding to a bulk Lorentz factor r re i = 1.1) relative to the counterstream. 
To better accommodate the parameters available for the calculations, it was necessary to make 
slight adjustments to the simulations we performed for comparison, leading to an increase in both 
the ion mass and, in the dilute case, the relative drift speed between the plasmas. The calculations 
also make predictions for modes at very large wavelengths in both spatial dimensions. We therefore 
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repeated the early stage of our simulation with the comparison parameters rrii/m e = 100 and 
r rel = 1.1, on an enlarged grid of 3840 A x 3840 A, and A sc = 30 A. 

To make a comparison with the analytical predictions, we extract the two-dimensional k spec- 
trum at times separated by only a short interval intended to capture the earliest linear growth of 
the instabilities, and compute the average growth rate, which is plotted in Figure [1] panels (c)-(e) 
and (h)-(j) for the dilute and symmetric cases, respectively. The agreement is satisfactory, both 
qualitatively and quantitatively, and the dominant modes are correctly rendered. As expected, the 
cold theoretical spectrum saturates at high k± while the hot version displays a local extremum for 
an oblique wave vector, as the kinetic pressure prevents the pinching of high-fc_i_ small filaments 



(|Silva et all 12002 ). 



Moreover, for high k±, the PIC spectrum describes wavelengths only a few cells long, small 
enough to be affected by the smoothing algorithm. This explains why these modes' growth is slower 
than expected. An electron skin depth of several hundred cells would almost certainly provide 
sufficient separation from the filtering length for the plots to agree better with the theoretical ones 
in this region, but such a simulation could not be large enough, or run long enough, to observe the 
later evolution without prohibitive computational expense. 

Modes with k± = are purely electrostatic, and produce no magnetic field. This is why their 
growth rate is much better rendered when measuring the E spectrum rather than the B spectrum. 
Conversely, modes with fen = are mostly electromagnetic with a very small phase velocity, which 
explains why their growth is only evident in the B spectrum. 



4. Results 

The qualitative behavior of the counterstreaming-plasma system in the absence of cosmic 
rays exhibits a dependence on the density ratio w particularly in the earliest stages of evolution, 
and at later times the impact of the magnetic field b becomes prominent. Although the details 
differ from one simulation to the next, the genera l behavior is similar to that seen in the three- 



dimensional simulations of iFrederiksen et al.1 (120041 ) . in which the collision of electron-ion plasmas 



is characterized by the formation of current channels first in the electrons, and later in the ions, 
which proceed to merge and grow. In our experiments, this growth of the channels and associated 
structure in the magnetic field eventually reaches a size comparable to the simulation domain, at 
which point the imposed periodic boundary conditions prevent further enlargement. However, the 
time necessary for this to occur — hundreds or thousands of uj~^ — may exceed the residence time of 
a given particle in the subshock. 

The benchmark behavior for our simulations of negligible-cosmic-ray energy density is plotted 
in Figure [2] (magnetic field absent), Figure [3] (weak magnetic field), and Figure 0] (strong magnetic 
field); the three columns of each plot correspond to the symmetric, intermediate, and dilute density 
ratios w. The role of the density ratio is prominent in the action of the two-stream instability on 
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the drift velocity of the counterstreaming electron populations (IMedvedev and Loebl Il999l ). For 
all considered values of magnetic field b, there is an abrupt deceleration of the electrons, both 
the stream and counterstream, when 20 < co pc t < 30. In the w = 1 symmetric case, this initial 
deceleration strips the electrons of nearly 90% of their relative drift, but in the dilute case, less than 
half the drift is removed. The electric and magnetic fields are amplified substantially at this point, 
but the reduced electron drift suppresses further immediate amplification. The ions are slower to 
respond, and their evolution differs qualitatively from that of the electrons: they form spatially 
alternating long-lived channels of current: as in lanes of vehicular traffic, ions moving one direction 
become spatially segregated from those moving the opposite direction, and this greatly lengthens 
the time for the ion drift velocities to converge, though considerable heating of the ions takes place 
prior to any significant systematic deceleration. 

The convergence of ion drift velocities reveals the primary effect of the magnetic field: the ion 
speeds remain well separated throughout the simulation lifetime in the absent or weak magnetic- 
field cases, but the strong magnetic field brings them together within roughly lO 4 ^" 1 for all density 
ratios. It may be that at least in this two-dimensional simulation, the transverse motion necessary 
for separating the ions into current channels is slightly inhibited by the guiding magnetic field, 
increasing the extent to which the counterstreaming ion populations are forced to interact. However, 
it is worth noting once more that the late-term behavior of the simulations is suspect on account 
of the structure size becoming comparable to the domain boundaries, an artificial upper limit. 



4.1. Behavior including cosmic rays 

For the configurations we considered, the presence of cosmic rays does not appear to result in 
any significant deviations from the behavior observed in their absence. When their energy density 
is given a value intended to represent conditions typical of the Galactic disk, no differences from 
the negligible-cosmic-ray configuration are observed. When cosmic rays are given an exaggerated 
abundance, subtle effects do appear in the simulation, but upon inspection they are disregarded 
for one of two reasons: they can be dismissed as numerical effects arising from the finite number of 
cosmic-ray particles employed in our simulation, or else they result in minor quantitative changes 
in the evolution that are of greatest prominence only when the effect of the periodic boundaries is 
already non-negligible. When the cosmic-ray weight is comparable to that of the plasma particles, 
their speed (approximately c) maximizes the amplitude of current-density fluctuations resulting 
from the statistically expected local departures from homogeneity. 

In order to verify that the observed differences resulted from our choice of representation and 
not from some underlying physics, we repeated an "abundant" simulation with a tenfold increase 
in computational particles representing the same physical cosmic-ray density. Figure [5] illustrates 
via the electromagnetic field amplitudes that the statistical noise levels arising at the earliest times 
saturate at a level vTO lower when cosmic rays are represented by 50 particles per cell instead 
of 5, bringing both the noise level and the detailed time evolution into better agreement with 
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the "negligible" case. A further significant increase in the computational particle count is too 
expensive for direct comparison with the simulations in Figure (5) Using a smaller computational 
grid, a simulation with 500 cosmic-ray particles per cell (leaving the other plasmas at their original 
20 per cell) illustrates the continuation of the trend observed with 50 per cell. Nevertheless, the 
remaining difference in non-noise behavior is already nearly imperceptible at just 50 computational 
particles per cell. This effect is paralleled in the other aspects of the system's evolution in which 
abundant cosmic rays appeared to result in minor differences, such as drift velocities and wave 
spectra. 



Discussion and conclusions 



Motivated by an interest in possible effects of cosmic rays on the physics governing the de- 
velopment of collisionless astrophysical shocks, we have performed several multidimensional PIC 
simulations of counterstreaming plasmas with various density ratios and magnetic-field strengths, 
both with and without a background population of energetic cosmic rays. This initially homo- 
geneous environment is intended to represent the interior of the subshock, or shock transition 
layer. Before cosmic rays are added to the picture, the system resembles the subject of numerous 
beam-plasma or interpenetrating-plasma studies, where the initial behavior of the system can be 
understood in terms of known instabilities. Most prominently, the counterstreaming electron pop- 
ulati ons are the first to interact , via symmetric or asymmetr ic two-stream instabilities, as seen in, 
e.g.. iMedvedev and Loebl (|1999l ) and iFrederiksen et al.l (|2004r). In particular, the thre e-dimensional 



simulations of unmagnetized electron-ion plasma collisions by IFrederiksen et al.l (|2004l ) demonstrate 
the formation of merging and growing current filaments in first the electrons and subsequently the 
ions, and the nuanced relationships among the various populations. 



One limitation of our approach is the reduced electron-ion mass ratio, m e /mj = 1/50. iBret and Dieckmann 



(|201Cj) have explored the effect of the mass ratio on the hierarchy of unstable modes in beam-plasma 
systems. By stating that a mass ratio different from 1/1836 should not change the nature of the 
most unstable mode, they were able to articulate a criterion defining the largest acceptable mass 
ratio. In the present case, the unstable linear spectrum depends very weakly on the mass ratio, 
when no cosmic rays are introduced. As stated in Section [3j finite-mass ions do not add any ex- 
tra unstable branches to the dispersion equation. As a result, the Bret-Dieckmann criterion is 
necessarily fulfilled in our configuration without cosmic rays. Turning now to the simulations in- 
cluding cosmic rays, we found no significant cosmic-ray effects with our current mass ratio within 
the cosmic-ray population. It is thus likely that a real mass ratio would result in an even weaker 
effect, leaving our conclusions unchanged. 

Even at just 50 times the electron mass, however, the behavior of the ions is markedly different. 
For the most part, the electrons produce and respond to turbulent small-scale electromagnetic fields 
that serve to mix their distribution functions; the electrons act in concert as a single population 
while the ions of the stream are still easily distinguished from the ions of the counterstream on 
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account of the filaments. Only when the filamentary structures have enlarged until relatively few 
repetitions are contained within the simulation plane do the ions begin to converge, and in some 
of our simulations that remains speculative on account of their finite duration. While it would 
be possible to extend the simulations ostensibly to observe the eventual convergence, this would 
be of limited value without simultaneously increasing the dimensions of the simulation domain to 
mitigate the effect of the periodic boundary conditions. An appreciable increase in size, however, 
may require us to abandon the simplicity of our present configuration by taking into account large- 
scale spatial variations, perhaps involving a clear distinction between upstream and downstream 
regions and an unstable charge-separation layer resulting from differences in the electrons' and ions' 
evolution. 

When we compare the system evolution in the presence of cosmic rays with that in their 
absence, we find that for the physical configurations studied, cosmic rays do not introduce a sta- 
tistically significant departure from the unperturbed results described in Section [5J This may be a 
consequence of the comparatively large mean free path and characteristic timescale for evolution of 
the cosmic-ray distribution: even with the amplification of electric and magnetic fields within the 
transition layer, cosmic rays of modest energy apparently do not couple to the dynamics of thermal 
electrons and ions in any appreciable way. We surmise that at least for unmagnetized or parallel 
subshocks, the impact of cosmic rays — even when their energy density is unusually large — on the 
instabilities mediating the subshock transition is negligible. 

Cosmic rays may still indirectly affect various pr operties of the shock , by modifying the up- 



stream environment from its quiescent characteristics (jStroman et al.ll2009l ): a shock will of neces- 
sity propagate differently through a turbulent, heated upstream medium — perhaps with a greatly 
amplified magnetic field — from the comparatively clean case of a uniform, cold interstellar medium 
in a gently fluctuating Galactic magnetic field. Irrespective of the cosmic-ray abundance, the rapid 
development of turbulence in the shock transition layer and the associated heating of the electrons 
in particular may provide an enlarged pool of candidates for injection into the standard diffusive 
shock acceleration mechanism. 

The combinations of parameters we explored did not yield any appreciable effects that could 
be attributed to the presence of cosmic rays. While we chose parameters intended to be relevant 
to nonrelativistic astrophysical shocks in environments where the presence of cosmic rays is sus- 
pected, it is nevertheless possible that in other, more exotic environments, cosmic rays may yet play 
some unforeseen role in the subshock microphysics. Future simulations in three dimensions and 
assigning additional degrees of freedom to the magnetic field and the cosmic-ray population may 
uncover effects that eluded the present analysis. However, a departure from the quasiparallel-shock 
configuration will also introduce effects not observed here, such as the development of magnetosonic 
waves, which can be more sensitive to the values of the mass ratio or oth er simulation parame ters; 



this in turn may affect the dissipation of excited unstable modes (see, e.g.. lHellinger et al.ll2007l ). In 
any case, the dearth of differences between even the grossly exaggerated cosmic-ray energy density 
and the system in which it was negligible provides a sense of reassurance that the physics of per- 
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haps a majority of astrophysical shock- forming instabilities can in principle be understood without 
invoking some direct microscopic interference by a spectator population of cosmic rays. 
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Fig. 1. — Measured or calculated growth rates (in units of electron plasma frequency) of unstable 
modes with wave vector ku , k± (scaled by X sc Av/c) for the dilute (left column) and symmetric (right 
column) stream configurations, with negligible-cosmic-ray density and the weak (b = 0.01) magnetic 
field. Within each column, the uppermost plot (a, f) is the instantaneous growth rate calculated 
in the zero-temperature limit. The next plot (b, g) includes a finite thermal spread of v t h_, e = 0.01c 
in the calculation. The remaining three rows are the average growth rate measured in the early 
stages of high-resolution PIC simulations, for the growth of perturbations in E x (c, h), E y (d, i), 
and B z (e, j). The growth-rate measurement is performed by calculating the fractional increase 
in the instantaneous two-dimensional Fourier power spectrum for the selected field component 
between times to and t\. For the dilute configuration, uj pc to ps 6 and u; pe ti ~ 18; for the symmetric 
configuration, u pc to ps 6 and w pc ii ps 12. By virtue of being an average over a finite interval, the 
growth rate obtained in this way is not a pure, instantaneous quantity like the top two rows of 
plots. The changing conditions in the stream and counterstream alter the instantaneous growth 
rate during the measurement interval, resulting in the minor differences that appear between the 
upper and lower parts of each column. 
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Fig. 2. — Selected results from the "negligible-absent" family of simulations, in which cosmic rays 
are negligible and an initial homogeneous magnetic field is absent. From left to right, the density 
ratio of the stream to the counterstream is w = 1, 0.3, and 0.1, corresponding to the "symmetric," 
"intermediate," and "dilute" cases, respectively. The bottom row shows the drift velocity of the 
stream electrons (solid line), stream ions (dot-dot-dot-dashed line), counterstream electrons (dot- 
dashed line), and counterstream ions (dashed line). A dotted line is included at zero velocity, 
representing the cosmic-ray ions and electrons. The second-to-bottom row displays the root-mean- 
square amplitudes of the electric and magnetic field parallel and perpendicular to the drift axis, 
where the magnetic field has been multiplied by c and all components are expressed in units of 
E u = w pc cm e /e. The second and third rows from the top illustrate the time evolution of the 
Fourier spectra of B z both parallel and perpendicular to the drift, with wavelengths measured 
in units of the electron skin depth A sc . A solid black line traces the evolution of the dominant 
values of \ x and X y . The top row plots the mean random kinetic energy (i.e., thermal energy) per 
particle, E rSLn = (7' — 1) mc 2 , where 7' is the particle Lorentz factor in the local plasma rest frame 
(measured over a square-shaped region of area 100A 2 ). The lines in the top row represent the same 
populations as in the bottom row. 
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Fig. 3. — Selected results from the "negligible-weak" family of simulations. The initial magnetic 
field is set such that the electron cyclotron frequency is f2 e = eB/m e = 0.01w pe . See the caption 
to Figure [2] for a detailed description of each plot. 
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Fig. 4. — Selected results from the "negligible-strong" family of simulations. The initial magnetic 
field is set such that the electron cyclotron frequency is f2 e = 0.1w pe . See the caption to Figure [2] 
for a detailed description of each plot. 
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Fig. 5. — Evolution of the electric and magnetic field components in the configuration b = 0.1, 
w = 0.3, as a function of the cosmic rays' computational representation. Panel (a) shows the 
case when cosmic rays are negligible. The "abundant" (u>cr7cr = 10) cosmic-ray distribution is 
represented by only five computational particles per cell in panel (c). In panel (e), the number 
of particles has been increased tenfold, and the initial field amplitudes decrease by approximately 
VTo, suggesting that statistical fluctuations in the cosmic-ray contribution to the local current 
density result in a high noise level in the field components. The column on the right reproduces 
the conditions of the left column on a grid whose dimensions have been reduced by 90% in each 
direction, allowing a full-length simulation with 500 cosmic-ray particles per cell. Panels (b) and (d) 
are comparable to panels (a) and (c), respectively, until uj pc t ~ 1000; the difference beyond this is a 
result of the smaller grid's periodic boundaries inhibiting the continued spatial growth of structures 
earlier. In panel (f), the initial noise level is reduced by a factor of \/l00 in all field components, 
corroborating our earlier supposition that the higher amplitudes are the effect of statistical noise. 



